import SimpleITK as sitk
import json
import numpy as np
from glob import glob
from pathlib import Path
import matplotlib.pyplot as plt

CASE_PATH = r'C:\Users\ang.a\OneDrive - Technion\Documents\MRI_Data\cases_for_rotem\99'
ANNOTATION_DATA = r'C:\Users\ang.a\OneDrive - Technion\Documents\MRI_Data\cases_for_rotem\99\data\24.11.2021__13_53..annotation.json'

# open and read the annotation data
with open(ANNOTATION_DATA, 'r') as f: 
    annotation_points = json.load(f)
    data = annotation_points['MPR points']

data = np.array(data)
data[-1] 

array([-11.43445065, -20.35332215, 101.2678986 ])

# transform the DICOM folder to a numpy array
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames(CASE_PATH)
reader.SetFileNames(dicom_names)
image = reader.Execute()
nda = sitk.GetArrayFromImage(image)

nda = np.flipud(nda)
nda = np.reshape(nda, nda.shape, order='F')

x_num = image.GetSize()[0]
xi = image.TransformIndexToPhysicalPoint([0,0,0])[0]
xf = image.TransformIndexToPhysicalPoint([x_num-1,x_num-1,0])[0]

y_num = image.GetSize()[1]
yi = image.TransformIndexToPhysicalPoint([0,0,0])[1]
yf = image.TransformIndexToPhysicalPoint([y_num-1,y_num-1,0])[1]

x_translate = ((np.abs(xi)+np.abs(xf))/2) - np.abs(image.GetOrigin()[0])
y_translate = ((np.abs(yi)+np.abs(yf))/2) - np.abs(image.GetOrigin()[1])

# get the slice indices that have annotations...
indices = [image.TransformPhysicalPointToIndex((0,0,z))[2] for z in np.unique(data[:,2])]
print(indices)

indices = list(range(indices[0], indices[-1]+1))
print(indices)

[47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]
[47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]

# convert the last coordinate to the slice index
data_copy = np.copy(data)

for i in data_copy:
    i[2] = image.TransformPhysicalPointToIndex((0,0,i[2]))[2]

i = 47
pt = (data[0])
size = image.GetSize()

plt.figure(figsize=(15,15))
plt.imshow(nda[i,:,:], extent=[xi,xf,yi,yf])
plt.plot(pt[0]+x_translate, pt[1]+y_translate, 'ro')

[<matplotlib.lines.Line2D at 0x232747448e0>]

# show all slices
fig, axs = plt.subplots(ncols=4, nrows=4, figsize=(50,50))

axs = np.ravel(axs)

for idx, i in enumerate(indices):
    axs[idx].set_title(f'SliceIdx: {i}')
    axs[idx].imshow(nda[i,:,:], extent=[xi,xf,yi,yf])

    for pt in data_copy:
        if pt[2] == i:
            axs[idx].plot(pt[0]+x_translate, pt[1]+y_translate, 'ro')
plt.tight_layout()